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Probabilistic Power Flow Computation via 
Low-Rank and Sparse Tensor Recovery 

Zheng Zhang, Hung Dinh Nguyen, Konstantin Turitsyn and Luca Daniel 


Abstract —This paper presents a tensor-recovery method to 
solve probabilistic power flow problems. Our approach generates 
a high-dimensional and sparse generalized polynomial-chaos 
expansion that provides useful statistical information. The result 
can also speed up other essential routines in power systems (e.g., 
stochastic planning, operations and controls). 

Instead of simulating a power flow equation at all quadrature 
points, our approach only simulates an extremely small subset 
of samples. We suggest a model to exploit the underlying low- 
rank and sparse structure of high-dimensional simulation data 
arrays, making our technique applicable to power systems with 
many random parameters. We also present a numerical method 
to solve the resulting nonlinear optimization problem. 

Our algorithm is implemented in MATLAB and is verified by 
several benchmarks in MATPOWER 5.1. Accurate results are 
obtained for power systems with up to 50 independent random 
parameters, with a speedup factor up to 9 x 10 20 . 

Index Terms —Power flow, power system, stochastic collocation, 
tensors, polynomial chaos, uncertainty, optimization. 

I. Introduction 

R EALISTIC power systems are affected by various uncer¬ 
tainties, such as the randomness of generations and loads, 
insufficient knowledge about network parameters, and noisy 
measurement m-m. Uncertainties may increase in future 
power systems, since many renewables highly depend on the 
uncertain weather conditions 0, @. These uncertainties must 
be considered in simulation, such that subsequent tasks can be 
completed in an efficient and robust way. 

This work investigates the probabilistic power flow prob¬ 
lem EP, which quantifies the uncertainties of bus voltages 
and line flows under uncertain loads, generations or network 
parameters. Currently, this problem is routinely solved in 
a number of decision-making procedures. Examples include 
transmission expansion and planning under long-term un¬ 
certainties in renewables penetration and regulation policies 
El, El- In operations, the operators assess the security of 
the system and calculate Available Transfer Capability using 
random scenario sampling ITTI where the ability to average the 
steady state solution over a large number of random scenarios 
is essential for secure power operations. 

Probabilistic power flow problems have been solved by 
Monte Carlo and many analytical methods (including multi¬ 
linearization El, the comulant method El, fuzzy load flow 
analysis El, and so forth). Recently, point estimation has 
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become a popular technique for probabilistic power flow 
analysis |4)-|[8]. This method assumes the solution being a 
summation of some univariate functions, then it computes the 
moments using a set of one-dimensional quadrature points. 

Stochastic spectral methods El have emerged as a promis¬ 
ing technique for the uncertainty analysis of many engineering 
problems including power systems El, El, El, EE 
They approximate the stochastic solution by a generalized 
polynomial-chaos expansion El- This representation can 
provide various statistical information (e.g., moments and 
probability density function); it can also accelerate many 
stochastic problems in power systems (e.g., stochastic unit 
commitment j9) and parameter inference 11221 ), whereas previ¬ 
ous approaches generally cannot. However, stochastic spectral 
methods may require lots of basis functions and simulation 
samples for problems with many uncertainties. In the uncer¬ 
tainty quantification community, some techniques based on 
compressed sensing (23), ||24| . proper generalized decomposi¬ 
tion (25], lt 26 l and tensor-train decomposition (20), (27], (28] 
have been developed for high-dimensional problems. 

This paper develops an alternative stochastic spectral 
method to solve probabilistic power flow problems with possi¬ 
bly high-dimensional random parameters. Our main contribu¬ 
tions are summarized as the following: i) We use tensors (29l 
(i.e., high-dimensional data arrays) to represent the huge set of 
data samples required in stochastic simulation. With a tensor 
format, we propose a low-rank and sparse tensor recovery 
scheme to generate a high-dimensional and sparse approxi¬ 
mation while using an extremely small subset of quadrature 
samples, ii) We present the detailed numerical implementa¬ 
tion of the tensor recovery method. Our algorithm relies on 
alternating minimization and the alternating direction method 
of multipliers (ADMM) (30). Although only locally optimal 
solutions are guaranteed, the developed solver performs well 
for many practical cases. We demonstrate the performance 
of the proposed technique with numerical simulations on 3 
benchmarks in MATPOWER 5.1 0. 

II. Problem formulation 
A. Probabilistic Power Flow Problem 

A steady-state power system with uncertainties can be 
described with parameterized power flow equations: 

n 

Pi{£) = J2 Vi(£)Vk(£) (Gikcos0ik{£) + f?ifcsin0ifc(£)) 
k =i 

n 

Q%(£)= Yj Vi(£)Vk(£) {Gik sm0ik(£) - Bik cos0ik(£)) 

k= 1 


where Pi, Qi, Vi, 9i are the active and reactive power, voltage 
magnitude and angle at load bus i, respectively; Gik and Ii lk 
are conductances and susceptances; da- = 0i~0k is the voltage 
angle difference between buses i and k. 

We employ random parameters £=[£i,-‘‘ , £d] £ R d to 
describe the uncertainties of load power consumptions that 
further influence bus voltages and angles. After computing 
Vi's and 0,'s, an out of interest y (e.g., the line flows) can be 
easily extracted. Obviously y also depends on £ and thus can 
be written as y = g (£). We assume that a deterministic solver 
is available to solve 0 given a sample of £. For simplicity, we 

assume that all elements of £ are mutually independent, then 

d 

their joint probability density function is p(£) = pk{£,k), 

fc =1 

where Pk{£,k) is the marginal probability density function of 
Moreover, the slack bus is assigned to compensate for the 
variations of loads and losses. 


B. Stochastic Collocation Method 

If the power flow problem is solvable, and y=g(£,) smoothly 
depends on £, then we can approximate y by a truncated 
generalized polynomial-chaos expansion Ell 

d 

y = g(£)~ <=„'!'<*(£), with <t>k, ak (^k). ( 2 ) 

| c*| <p k= 1 


The multivariate polynomial basis v k a (£) is indexed by 
«=[«!,••• ,ad]sN d , with the total polynomial degree 

d 

l«l=£ I a k | < p. The total number of basis functions is 

k= 1 


= (p + d)\ 
p\d\ 


(3) 


As shown in Appendix [A] the degree-a j. univariate polyno¬ 
mials {4>k,a k {£,k )} P ak —o are orthonormal to each other. There¬ 
fore, the multivariate basis functions are also orthonormal, and 
c a can be computed with projection 


c a = f ^ a {€)g{€)p{€)d€- (4) 

R d 


This integral can be evaluated by a proper quadrature rule 
which requires computing g(£) at a set of samples. 



Fig. 1. Demonstration of vectors (left), matrices (middle) and tensors (right). 


III. A Tensor-Recovery Approach 

This section presents our tensor-recovery method to solve 
high-dimensional probabilistic power flow problems. 


A. Tensor Representations of 0 

As a generalization of vectors and matrices, a tensor A. £ 
I" 1 ’ x"-xmd represents a high-dimensional data array (29l . 
The number of dimensions, d, is called the mode of a 
tensor; mk is the size of the k- th dimension. Given index 
i = (*!,■■• ,id) (with integer ik £ [1,to/c]), we can specify 
one element *4.(i). Fig.Q]shows a 1-mode tensor (i.e., vector), 
a 2-mode tensor (i.e., matrix) and a 3-mode tensor. 

First, we define a d-mode tensor Q £ Ir x '" xm 

00) = s(£v ■■,#)• ( 6 ) 

Next, for every £& and its degree-a*,- polynomial (f>ka k {tik), 
(k) 

we define a vector w ak £ R m with its ik -th element being 

w a* {ik) = (t>k, ak (Ck ) w k ■ ( ? ) 

For every index vector a, we further construct a d-mode rank- 

1 tensor W a £ M mx - xm : 

d 

= wW o • • • o W<J W a (i) = Yl w a k {ik)- (8) 

fc=l 

Here o denotes an outer product. As a result, the right-hand 
side of 0 is the inner product of Q and W„: 


C. Integration Rules and Curse of Dimensionality 

Among different quadrature rules Q3-0D, this work 
considers computing c a by a tensor-rule Gauss quadrature 
method. First, use Gauss quadrature f35l (in Appendix 10 to 
decide m quadrature samples and weights {Ck ; w k }, , f° r 

£k- Next, we compute c a by a tensor rule 

mm d 

c “~ ’Q) ■ (5) 

i\=l id =1 k —1 

This method requires simulating the power flow equa¬ 
tion m d times, and obviously it only works well for low¬ 
dimensional problems (e.g., when d is below 5 or 6). Sparse 
grid has been applied to simulate power systems m, which 
can compute 0 with about 2 P K samples for high-dimensional 
cases lfl9l . In this paper, we aim to use only < K samples 
from a tensor rule to compute 0. 


Ca « (g, w„> = ^ • • • ]T Q{ i)w Q (i). 

*1 = 1 id = 1 


(9) 


In summary, in order to obtain the generalized polynomial- 
chaos approximation 0 we need to compute: 1) tensor Q ; 2) 
tensor W„ for each a. satisfying |a| < p. Since each W a is 
the outer product of d vectors and many of them can be reused, 
computing W a ’s is trivial. However, directly computing Q 
is almost impossible, since the power flow equation must be 
simulated m d times. 


B. Low-Rank and Sparse Tensor-Recovery 

Instead of computing Q directly, we approximate Q by 
tensor recovery. The key idea is described below. 





1) Sub-Sampling: We randomly compute a small portion of 
elements in Q , then seek for a tensor Q to approximate Q. Let 
I = {i|l < ffc < m} include the indices for all elements in Q. 
The size of X, \X\, is m d . We choose a subset f 1 C X (with 
|f2| -C \X\) that includes a small number of indices randomly 
selected from X, and compute Q{ i) = g(Ci i ■ ■ ■ , C/) f° r an y 
i £ Cl. Then, we look for a tensor Q such that it matches Q 
at all elements specified by Cl, i.e., 

||Pn (q - g) III = 0. (10) 

Here Pn is a linear operator for tensors: 


B = P n (A) B( i) 


«4.(i), if i £ Cl 
0, otherwise. 


The Frobenius-norm of a general tensor is defined as 


II^IIIf = ViAA). (12) 


An infinite number of tensors exist that satisfies the require¬ 
ment ([Tol l but significantly differs from Q. Therefore, some 
constraints can be added to regularize this problem. 

2) Constraint 1- Sparsity: Let vector c = [■ ■ ■ , c a , ■ ■ ■ ] £ 
R K includes all coefficients in the generalized polynomial- 
chaos approximation. In high-dimensional cases, c is generally 
very sparse - most of its elements are close to zero. Using 1/ - 
norm as a measure of sparsity l36l , we have 


CK<p 


E 

a<p 




is small 


(13) 


3) Constraint 2-Low Tensor Rank: In many cases, Q has 
a low tensor rank and can be well approximated by the 
summation of a few rank-1 tensors. Therefore, we assume that 
the solution Q has a rank-/’ decomposition: 

S = T (U (1) , • • • , U (d) ) := uf } o • • • o (14) 

3=1 


where ir: fc) £R mxl is the j-th column of matrix Ul fc l £ R mxr . 
Therefore, we may use d matrices to represent the whole tensor 
instead of computing and storing all elements of Q. 

4) Final Tensor-Recovery Model: We describe the low-rank 
and sparse tensor-recovery model as follows: 

Given <?(i) for every i £ Cl, solve 


min /(uw,- ,u (d) ) 

•fuW gR mxr }' i _ 

= i||Pa (T(U«,- ,U( d ))-a)||^ 
+A E |(T(UW,--- ,UW),W Q )|. 

M<p 


(15) 


Here A > 0 is a regularization parameter. 


C. Summary of Main Steps 

We summarize the main steps of our approach as below. 

1) Simulation Step: Randomly generate a small subset 
Cl C X such that |fl| < K <C m d . For every index 
i = [*!,••■ ,id\ £ Cl, simulate the power flow equation <[TJ 
once to obtain a deterministic value £?(i) = p(^ 1 , • • • ,Cd)- 


Algorithm 1 Alternating Minimization for Solving (IT5T >. 

l: Initialize: £ R mx? ' for k = 1, • • • d; 

2 : for l = 0, 1, • • • 

3: for k = 1, • • • , d do 

4: solve dm by Alg.[2]to obtain uT fc ) ,z+1 ; 

5: end for 

6: break if converged; 

7: end for 

8: return = Ul fe l ,!+1 for k = L • • ■ ,d. 


2) Optimization Step: Solve d to obtain matrices 
U*- 1 ), • • • , that represent tensor Q in m. 

3) Model Generation: Replace Q by Q, and calculate c a ’s 
according to (O. With low-rank tensor factors, the computation 
can be simplified to 


<?> W «) = E 


3 =1 


a / 

n((“fy 

,fe=i v 


,(*) 


(16) 


which involves only cheap vector inner products. 

Since we can approximate Q by using only a small number 
of simulation samples, our method can be applied to many 
high-dimensional problems. 


IV. Optimization Solver 
This section describes how to solve (IT5l >. 


A. Outer Loop: Alternating Minimization 

1) Algorithm Flow: Starting from an initial guess 

{U(fc), 0 }d = 

.j, we perform the following iterations: at iteration 
( + 1 we use as an initial guess and obtain updated 

tensor factors {Ul fc l ,i+1 }^_ 1 by alternating minimization. 
Each iteration consists of d steps, and at the fc-th step, 
is obtained by solving 




arg min / 



(fc-l),Z+l 


,X,U (fe+1) ’ i 



Since all factors expect are fixed, d becomes a 

convex optimization problem, and its global minimum can 
be computed by the solver in Section IIV-BI The alternating 
minimization method ensures that the cost function decreases 
monotonically to a local minimal. The pseudo codes are 
summarized in Alg.[T] which terminates when the convergence 
criteria in Appendix ICl is satisfied. 

2) Prediction Error: An interesting question is: how ac¬ 
curate is Q compared with the exact tensor £?? Our tensor 
recovery formulation enforces consistency between Q and Q 
at the indices specified by O. We hope that Q also has a good 
predictive behavior - Q(i) is also close to i) for i ^ Cl. 
In order to measure the predictive property of our results, we 
define a heuristic prediction error 


c p r 


\ 


E (c(i)-e(i 

ieO' v 


e mwm 

iGfi' 


with W\ = ]^[ (18) 


k =1 

















Algorithm 2 ADMM for Solving ([UT >. 

1 : Initialize: form A,F and b according to Appendix iDl 
specify initial guess x°, u° and z°; 

2 : for j = 0,1, • • • do 

3: compute x J+1 , z J+1 and u - ?+1 according to (f20t ; 

4: break if ||Fx J+ 1 —z J+1 || < ei & ||F T (z J+ 1 —z J )|| < 62 ; 

5: end for 

6 : return U ^’ i+1 = reshape(x J+1 , [to, r]) . 



Here 11' £l is a small-size index set such that SI' D O = 0. 
Obviously, Q has good predictive behavior if e pr is small. 
Estimating e pr requires simulating the power flow equation at 
some extra quadrature samples. However, a small-size Cl' can 
provide a good heuristic estimation. 


B. Inner Loop: Numerical Solver for ( I / 71 ) 

Following the procedures in Appendix [D] we rewrite Prob¬ 
lem (fT71l as the generalized LASSO problem: 


^Tj(fc),(+i 


= arg mm ■ 


I Ax — b|| 2 + A|Fx| 


(19) 


where A e Rl n l xmr , F e R Kxmr and b <E K |n|xl , and 
x = vec(X) € M mrxl is the vectorization of X (i.e., x(jm — 
m + i) = X(i,j) for any integer 1 < i < m and 1 < j < r ). 
Note that |S7| is the number of simulations samples in tensor 
recovery, and K is the total number of basis functions. 

We solve ( fl~9b by the alternating direction method of mul¬ 
tipliers (ADMM) ll30l . Problem ( 1 1 91 can be rewritten as 
1 9 

min — ||Ax — b|| 2 + A|z| s.t. Fx — z = 0. 

X,Z 2 


By introducing an auxiliary variable u and starting with initial 
guesses x°, u° = z° = Fx°, the following iterations are 
performed to update x and z: 

x fc+1 = (A t A + sF t F ) _1 (A T b + sF T (z fe - u fc )) 

z k+1 = shrink A/s (Fx fc+1 + z k + u k ) (20) 

u k+i = u fc + Fx fc+i _ z k+i' 


Here s > 0 is an augmented lagrangian parameter, and the 
soft thresholding operator is defined as 

{ a — A/s, if 0 > A/s 
0, if |a| < A/s 
a + A/s, if 0 < —A/s. 

The pseudo codes for solving (IUI i are given in Alg. [ 2 ] 


C. Limitations 

Firstly, the cost function of (IUI i is non-convex, and it is 
non-trivial to compute its global minimum with theoretical 
guarantees. Although researchers and engineers are very often 
satisfied with a local minimal, the obtained result may not 
be good enough for some cases. Secondly, in this work the 
parameters A [the regularization parameter in ( fUl il and r 
[the tensor rank in (fUT il are set based on some heuristic 
experiences. This treatment is definitely not optimal and does 
not guarantee high accuracy for all cases. 


Fig. 2. The 6-bus system. 




Fig. 3. Left: all quadrature samples; right: samples used in tensor recovery. 


V. Simulations 

This section reports the simulation results for several test 
cases from MATPOWER 5.1 ED. All codes are implemented 
in MATLAB. We find that a 2nd- or 3rd-order generalized 
polynomial-chaos expansion can provide good accuracy for 
many cases, therefore we set p = 2 (or 3) in <f2]) and m = 3 
(or 4) in Equation (|9} . 

A. 6-Bus Case (with 3 Random Parameters) 

The case6ww example in MATPOWER 5.1 (c.f. Fig. 0 is 
used as a demonstrative example. We use 3 random parameters 
to describe the uncertain active powers at the load buses 4 
to 6 . We aim to obtain a 3rd-order generalized polynomial- 
chaos expansion for the real power injected from Bus 2 to 
Bus 4, leading to 20 basis functions in total. Applying a 4- 
point Gauss-quadrature rule to perform numerical integration 
for each dimension, we generate 64 quadrature points in total. 

In order to compute the generalized polynomial-chaos ex¬ 
pansion, only 18 quadrature points (as shown in Fig. 0 are 
randomly sub-selected. The simulation results at these selected 
samples are used to perform tensor recovery. For this case, we 
find that setting the tensor rank r = 3 and the regularization 
parameter A = 0.25 is a good choice. Starting from a randomly 
generated rank-3 tensor, our algorithm converged after 25 
iterations as shown in Fig. Q] The obtained low-rank tensor 
approximation has an estimated prediction error of 0 . 2 %. 
With the obtained tensor approximation, the coefficients for 
all generalized polynomial-chaos basis functions are easily 
calculated based on (f]~6l >. The coefficients for a = 0 is 31.83, 
which is the mean value of the output. All other coefficients 
are plotted in Fig. 0 where a sparsity pattern is observed. 

Next we validate our results by Monte Carlo. Here we 
use 5000 samples in Monte Carlo simulation and treat its 
result as a golden reference solution. As shown in Table [j] the 



















Fig. 4. Convergence for the 6-bus example. The left figure shows the relative update of tensor factors (i.e., extensor); the middle figure shows the update 
of the generalized-polynomial-chaos coefficients (i.e., e; ig pc); the right figure shows the value of the objective function (i.e., //). 
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Fig. 5. Coefficients of {'J'o;(4)}^ a |—j tor the 6-bus example. 


TABLE I 

Moments of the 6-bus example. 



Tensor Recovery 

Monte Carlo 

samples 

18 

5000 

Mean 

31.83 

31.87 

stand, dev. 

0.0439 

0.0448 


mean value and standard deviation from our tensor recovery 
approach is very close to that from Monte Carlo. 

Complexity Reduction. Since we use 18 samples out of 
64 quadrature points, the reduction ratio for this problem is 
3.6. Note that the number of samples in tensor recovery is less 
than the number of basis functions (i.e., 20). 

B. 30 -Bus Case (with 24 Random Parameters) 

Next we consider the case30 example in MATPOWER 
5.1, with the active powers of 24 load buses modeled by 
Gaussian random variables. We apply a 2nd-order generalized 
polynomial-chaos expansion for the real power from bus 15 to 
bus 23, requiring totally 325 basis functions. For each param¬ 
eter, 3 quadrature points are used, leading to 3 24 ~ 2.8 x 10 11 
samples in total. Obviously, it is prohibitively expensive to 
simulate the power system at all quadrature points. 

In our tensor recovery scheme, we randomly pick 280 
quadrature points from the full tensor-rule quadrature samples 
and approximate Q by a rank-4 tensor. Setting A = 0.3 and 
starting with a random initial guess, our algorithm converges 
nicely’ after 26 iterations which are similar to Fig. ED With 50 
newly sub-sampled quadrature points as the testing samples. 
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Fig. 6. Coefficients of {'l'ct(£)} 2 Q; |_ 1 for the 30-bus example. 
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Fig. 7. Flistograms of the simulated output for the 30-bus example . 


the estimated prediction error is 0.55%. Although this example 
has many random parameters, its generalized polynomial- 
chaos expansion is very sparse, as shown in Fig. [6] 

In order to check the accuracy, we perform Monte Carlo 
simulation using 5000 random samples. Table [II] compares 
the mean values and standard deviations from both ap¬ 
proaches, and they are very close. An advantage of generalized 
polynomial-chaos expansion is that one can easily evaluate the 
expression with many samples to get a density function or 
histogram. Such information cannot be easily obtained by a 
point-estimation method. The histogram from our method is 
close to that from Monte Carlo (c.f. Fig. [7]). 

Complexity Reduction. Since we use 280 samples out of 
3 24 quadrature points, the reduction ratio for this example is 
10 9 . The number of samples in tensor recovery is also smaller 
than the number of basis functions (i.e., 325). 
































10 3 





Fig. 8. Convergence for the 57-bus example. The left figure shows the relative update of tensor factors (i.e., e;,tensor); the middle figure shows the update 
of the generalized-polynomial-chaos coefficients (i.e., ej !g pc); the right figure shows the value of the objective function (i.e., ft). 



Fig. 9. Coefficients of tor the 57-bus example. 


TABLE II 

Moments of the 30-bus example. 



Tensor Recovery 

Monte Carlo 

samples 

280 

5000 

Mean 

-10.23 

-10.22 

stand, dev. 

0.048 

0.049 


C. 57-Buse Case (with 50 Random Parameters) 

Finally we consider the case57 example in MATPOWER 
5.1, with 50 Gaussian random variables describing the active 
powers at load buses. With a 2nd-order polynomial-chaos 
expansion, we aim to approximate the real power injected 
from Bus 19 to Bus 20 with 1326 basis functions. Using 
3 Gauss-quadrature points for each parameter, a tensor-rule 
quadrature method requires 3 50 > 7 x 10 23 samples in total. 
It is impossible to store the samples on a personal computer, 
let alone simulating the power flow equation at all samples. 

Our tensor recovery scheme randomly sub-selects 800 
samples to perform power flow simulations. Starting with a 
random initial guess, we approximate the full tensor Q by 
a rank-5 tensor, with an estimated prediction error of 1%. 
Fig. [8] shows the convergence of our solver. Fig. [9] plots the 
coefficients for all non-constant basis functions. Clearly, the 
result is extremely sparse for this high-dimensional example. 

In order to get a full picture about the statistical behavior of 
the output, we evaluate the computed generalized polynomial- 
chaos expansion with 5000 random samples and plot its proba¬ 
bility density function. As shown in Fig. [To] the result is close 
to that from Monte Carlo simulation on the original power flow 


Fig. 10. Computed probability density functions for the 57-bus example. 


TABLE III 

Moments of the 57-bus example. 



Tensor Recovery 

Monte Carlo 

samples 

800 

5000 

Mean 

0.721 

0.724 

stand, dev. 

0.0596 

0.0609 


equations. The mean values and standard deviations from both 
approaches are very close (c.f. Table HIiV 

Complexity Reduction. Since we use 800 samples out of 
3 50 quadrature points, the reduction ratio for this example is 
about 9 x 10 2 °. The number of samples in tensor recovery is 
again smaller than the number of basis functions (i.e., 1326). 

VI. Conclusions and Future Work 

This paper has presented a probabilistic power flow simula¬ 
tion algorithm based on tensors and stochastic collocation. In 
order to break the curse of dimensionality, we have developed 
a high-dimensional method that exploits the low-rank and 
sparse property of tensors. This tensor framework has com¬ 
pleted the huge-data-size simulation task with an extremely 
small-size simulation data set. We have further developed a 
numerical solver for the tensor recovery problem and tested 
it on three power flow benchmarks. This algorithm has suc¬ 
cessfully generated high-dimensional and sparse generalized 
polynomial-chaos expansions for the solutions. Good accuracy 
(measured by prediction errors and comparison against Monte 
Carlo) as well as significant computational cost reduction (with 
up to 9 x 10 2 ° times) have been observed in this work. 







































To the best of our knowledge, this is the first work studying 
stochastic power systems from a tensor perspective. While 
several limitations exist in the current work, the authors believe 
that many problems are worth investigation in this direction. 
Firstly, it is worth investing global non-convex optimization 
to solve m °r trying to convexify the formulation. Secondly, 
some future work about the optimal choice of A and r may 
help improve the robustness of our framework. Finally, it is 
worth developing better sampling schemes to improve the 
performance of tensor recovery. 

One particularly attractive application of this technique is 
the construction of probabilistic static equivalents of a sub¬ 
networks. These equivalents can be used both for distribution 
grid models with high penetration of intermittent renewables 
and for probabilistic modeling of random disturbances in 
neighbor areas in multi-area power systems. 

A second natural application is the stochastic contingency 
analysis over a range of operating conditions. Currently de¬ 
terministic contingency analysis is considered in some basic 
cases with given operating condition. However, as the systems 
change and are subject to uncertainties, one need to take 
multiple cases or scenarios into account. The fast averaging 
technique developed in this paper can effectively alleviate the 
heavy computation in such situations. 


Appendix A 

Orthonormal Polynomials 
Consider a single random parameter ^ el with a prob¬ 
ability density function pk{£,k), one can construct a set of 
polynomial functions subject to the orthonormal condition: 

4>k,a{^k)(/}k,l3(£,k)Pk{(k)d£,k = $aj3 
R 

where 5 a> p is a Delta function, integer a is the highest 
degree of 4> k ,a(t;k)- Such polynomials can be constructed 
as follows [37]|. First, one constructs orthogonal polynomials 
{7 r fe,a(^fe)}a=o w itb an leading coefficient 1 recursively 

7I~fc,a+l ) — (£fc la) ^T,ck(Qc) K , a'Kk,ac—l(.^k) 

for a = 0, 1, with initial conditions TTk.-i(£, k ) = 

0, 7^0 (£fc) = 1 and = 1- For a > 0, the recurrence 
parameters are defined as 

E (W Q fe)) E(W q+ i(£ 

la = -7-x-, Na+1 = -7--- 

E ( 7r fe,a(^)) E (& 7r fc,a(& 

Here E denotes the operator that calculates expectation. Sec¬ 
ond, one can obtain {</>&,a (£&)}«=o by normalization: 

<? h,a(€k ) = — , for a = 0,1, • ■ ■ ,p. 

y/K 0 Ki ■■■ n a 



Appendix B 

Gauss Quadrature Rule |[35l 
Given ^ gl with a density function pk(£,k) and a smooth 
function q(£ k ), Gauss quadrature evaluates the integral 


/ 


q(tk)Pk(tk)d£k ~ diCkWk 

ik =1 


with an error decreasing exponentially as to increases. An 
exact result is obtained if q(£, k ) I s a polynomial function 
of degree < 2 to — 1. One can obtain {{C k > wl k )}i + =i by 
reusing the recurrence parameters in (EH to form a symmetric 
tridiagonal matrix J 6 r(p+ 1 ) x (p+ 1 ) : 


J (i,j) = < 


li—i i if i=j 
yTcI, if i = j + 1 
^/TcJ, ifi=j-l 
0, otherwise 


for 1 < i, j < p + 1. 


Let J = QLQ 7 be an eigenvalue decomposition and Q a 
unitary matrix, then ^ = Y\(i k ,i k ) and w % £ = (Q(l,ifc)) • 


Appendix C 

Error Control in Alg.Q] 

With tensor factors obtained after l iterations 

of the outer loops of Alg. [Q we define 

fi := / [updated cost func. of (fTSll I 

Gi := T • • • , (approximated tensor) 

cL := (Gi, Wo) [updated coefficient for 'F a (a)[ 

and let c l = [• • • , c l a , • • • ] £ for all |ck| < p. Then, we 
define the following quantities for error control: 

• Relative update of the tensor factors: 


CZ,tensor = . 5Z H U(fc) ’' ~ UW-*- 1 |||/^ ||U<*).*-1|| 
\ fc= 1 


- 1112 ^ 


*:=i 


Relative update of c = 


Q,gPC = 


= C - c 


Z —11 


„Z-li 


• Relative update of the cost function: 

Q,cost = l/z fl — 11/|/z — 11 • 

The solution is regarded as a local minimal if 

ez,tensor, £z, g pc and e ZjCOSt are small enough. 


Appendix D 

Assembling The Matrices and Vector in (fl9l i 

Consider the tensor factors U( fc_1 b i+1 , X, 

U( fe+1 ) ,i , ■■■, ffi ( |T7l ). We denote the (i,j) element 

of Tj( fc (or X) by u[ k d' 1 (or Xi j), and its j-th column by 

Uj k )J (or Xj). Then, the cost function in ( 1 1 7| i is 


/(- 


■■■ ,U' 


(k-l),l + l 


,X,U (fe+1) ’ ! ,.--) 


ien \j= 1 / |a|<p j= 1 

where the scalars // i ( and v a j are computed as follows: 


k -1 


mm= n <>r n 

k f =k+i 


( k'),l 

u i k , ,j , 


k '=1 
k -1 




k '=1 


k'=k+l 











Since each row (or element) of A (or b) corresponds to 
an index i £ SI, and each row of F corresponds to a basis 
function '!•'„(£;), in this appendix we use i as the row index 
(or element index) of A (or b) and a. as the row index of F. 
Now we specify the elements of A, b and F of ( fl9] >. 

• For every i £ SI, b(i) = <7(i). 

• Since Xi k j is the ( j — 1 )m + ik- th element of x = 
vec(X), for every i £ SI we have 


A(i, (j - 1 )m + i k ) 


Mi ,j, for j = 1, • • • ,r 
0, otherwise. 


■ Since x y includes the elements of x ranging from index 
(j — l)m + 1 to jm, given an index vector a the 
corresponding row of F can be specified as 

F (ajm-m + i k ) = ^jW^ik) = v a ,j4>k, ak (CkWk 


for all integers j £ [l,r] and ik £ [1 , m]. 


References 

[1] A. J. Conejo, M. Carrion, and J. M. Morales, Decision making under 
uncertainty in electricity markets. Springer, 2010, vol. 1. 

[2] B. Borkowska, “Probabilistic load flow,” IEEE Trans. Power Apparatus 
and Syst., vol. 93, no. 3, pp. 752-759, 1974. 

[3] E. Haesen, J. Driesen, and R. Belmans, “Stochastic, computational and 
convergence aspects of distribution power flow algorithms,” in Proc. 
IEEE Power Tech. , 2007, pp. 1447-1452. 

[4] C. S. Saunders, “Point estimate method addressing correlated wind 
power for probabilistic optimal power flow,” IEEE Trans. Power Syst., 
vol. 29, no. 3, pp. 1045-1054, 2014. 

[5] J. M. Morales and J. Perez-Ruiz, “Point estimate schemes to solve the 
probabilistic power flow,” IEEE Trans. Power Syst., vol. 22, no. 4, pp. 
1594-1601, 2007. 

[6] J. M. Morales, L. Baringo, A. J. Conejo, and R. Mrnguez, “Probabilistic 
power flow with correlated wind sources,” IET Generation, Transmission 
& Distribution, vol. 4, no. 5, pp. 641-651, 2010. 

[7] G. Verbic and C. Canizares, “Probabilistic optimal power flow in 
electricity markets based on a two-point estimation method,” IEEE 
Trans. Power Syst., vol. 21, no. 4, pp. 1883-1893, 2006. 

[8] C. Su, “Probabilistic load-flow computation using point estimate 
method,” IEEE Trans. Power Syst., vol. 20, no. 4, pp. 1843-1851, 2005. 

[9] E. M. Constantinescu, V. M. Zavala, M. Rocklin, S. Lee, and M. An- 
itescu, “A computational framework for uncertainty quantification and 
stochastic optimization in unit commitment with wind power genera¬ 
tion,” IEEE Trans. Power Syst., vol. 26, no. 1, pp. 431—441, Feb. 2011. 

[10] R. Allan and A. Leite da Silva, “Probabilistic load flow using multi- 
liearizations,” IET Gen. Transm. Distrib., vol. 128, no. 5, pp. 280-287, 
1981. 

[11] R. Allan and M. Al-Shakarchi, “Probabilistic techniques in AC load flow 
analysis,” Proc. 1EE, vol. 124, no. 2, pp. 154-160, Feb 1977. 

[12] V. Miranda, M. Matos, and J. Saraiva, “Fuzzy load flow: new algorithms 
incorporating uncertain generation and load representation,” in Proc. 
Power Syst. Comp. Conf., 1990, pp. 621-627. 

[13] G. Lin, M. Elizondo, S. Lu, and X. Wan, “Uncertainty quantification in 
dynamic simulations of large-scale power system models using the high- 
order probabilistic collocation method on sparse grids,” Int. J. Uncert. 
Quant., vol. 4, no. 3, pp. 185-204, Feb. 2014. 

[14] J. Hockenberry and B. Lesieutre, “Evaluation of uncertainty in dy¬ 
namic simulation of power system models: the probabilistic collocation 
method,” IEEE Trans. Power Syst., vol. 19, no. 3, pp. 242-272, 2004. 

[15] M. Milligan, P. Donohoo, and M. OMalley, “Stochastic methods for 
planning and operating power system with large amounts of wind and 
solar power,” in Proc. Int. Workshop Large-Scale Integr. Wind Power 
into Power Syst., 2012. 

[16] B. Gorenstin, N. Campodonico, J. da Costa, and M. Pereira, “Power 
system expansion planning under uncertainty,” Power Systems, IEEE 
Transactions on, vol. 8, no. 1, pp. 129-136, Feb 1993. 

[17] J. Stahlhut, F. Gao, K. Hedman, B. Westendorf, G. Heydt, P. Sauer, 
and G. Sheble, “Uncertain power flows and transmission expansion 
planning,” in Power Symposium, 2005. Proceedings of the 37th Annual 
North American. IEEE, 2005, pp. 489^496. 


[18] D. Xiu, “Fast numerical methods for stochastic computations: A review,” 
Commun. Comput. Phys., vol. 5, no. 2-4, pp. 242-272, Feb. 2009. 

[19] Z. Zhang, T. A. El-Moselhy, I. M. Elfadel, and L. Daniel, “Stochastic 
testing method for transistor-level uncertainty quantification based on 
generalized polynomial chaos,” IEEE Trans. Computer-Aided Design 
Integr. Circuits Syst, vol. 32, no. 10, pp. 1533-1545, Oct 2013. 

[20] Z. Zhang, I. Osledets, X. Yang, G. E. Kamiadakis, and L. Daniel, 
“Enabling high-dimensional hierarchical uncertainty quantification by 
ANOVA and tensor-train decomposition,” IEEE Trans. CAD Integr. 
Circuits Syst., vol. 34, no. 1, p. 63 76, Jan. 2015. 

[21] D. Xiu and G. E. Kamiadakis, “The Wiener-Askey polynomial chaos for 
stochastic differential equations,” SIAM J. Scientific Computing, vol. 24, 
no. 2, pp. 619-644, Feb 2002. 

[22] Z. Zhang, C. Petra, N. Petra, and E. Constantinescu, “Fast state esti¬ 
mation of power systems with polynomial chaos,” Aug. 2015, technical 
Report, Argonne National Laboratory. 

[23] X. Yang and G. E. Kamiadakis, “Reweighted l\ minimization method 
for sothcastic elliptic differential equations,” J. Comp. Phys., vol. 248, 
no. 1, pp. 87-108, Sept. 2013. 

[24] J. Peng, J. Hampton, and A. Doostan, “A weighted Zi minimization 
approach for sparse polynomial chaos expansion,” J. Comp. Phys., vol. 
267, no. 1, pp. 92-111, Jun. 2014. 

[25] A. Nouy, “Proper generalized decomposition and separated representa¬ 
tions for the numerical solution of high dimensional stochastic prob¬ 
lems,” Arch. Comp. Meth. Eng., vol. 27, no. 4, pp. 403-^134, Dec 2010. 

[26] F. Chinesta, P. Ladeveze, and E. Cueto, “A short review on model order 
reduction based on proper generalized decomposition,” Arch. Comput. 
Meth. Eng., vol. 18, no. 4, pp. 395^-04, 2011. 

[27] B. N. Khoromskij and C. Schwab, “Tensor-structured Galerkin approxi¬ 
mation of parametric and stochastic elliptic PDEs,” SIAM J. Sci. Comput, 
vol. 33, no. 1, pp. 364-385, Oct 2011. 

[28] D. Bigoni, A. P. Engsig-Karup, and Y. M. Marzouk, “Spectral tensor- 
train decomposition,” arXiv preprint, arXiv:1405.5713vl, May 2014. 

[29] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” 
SIAM Review, vol. 51, no. 3, pp. 455-500, Aug. 2009. 

[30] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed 
optimization and statistical learning via the alternating direction method 
of multipliers,” Foundations and Trends in Machine Learning, vol. 3, 
no. 1, pp. 1 -122, 2010. 

[31] R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas, “MAN¬ 
POWER: Steady-state operations, planning and analysis tools for power 
system research and education,” IEEE Trans. Power Syst., vol. 26, no. 1, 
pp. 12-16, Jun. 2011. 

[32] T. Gerstner and M. Griebel, “Numerical integration using sparse grids,” 
Numerical Algorithms, vol. 18, pp. 209-232, Mar. 1998. 

[33] S. Weinzierl, “Introduction to Monte Carlo methods,” NIKHEF, Theory 
Group, Amsterdam, The Netherlands, Tech. Rep., 2000. 

[34] W. J. Morokoff and R. E. Caflisch, “Quasi-Monte Carlo integration,” J. 
Comput. Phys., vol. 122, no. 2, pp. 218-230, Dec 1995. 

[35] G. H. Golub and J. H. Welsch, “Calculation of Gauss quadrature rules,” 
Math. Comp., vol. 23, pp. 221-230, 1969. 

[36] D. L. Donoho, “Compressed sensing,” IEEE Trans. Informa. Theory, 
vol. 52, no. 4, pp. 578 -594, April 2006. 

[37] W. Gautschi, “On generating orthogonal polynomials,” SIAM J. Sci. Stat. 
Comput., vol. 3, no. 3, pp. 289-317, Sept. 1982. 


